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Abstract Previous work found that different Japanese lineages of salamanders had quite different levels of species and 
genetic diversity. Lineages vary from having one to several species and the extent of genetic variation among lineages 
differs substantially. Most speciose, genus Hynobius contains 18 species and several potential cryptic species. We 
explore genetic diversity in this genus by combining comprehensive sampling and mitochondrial DNA sequences. Based 
on this and previous analyses of salamanders, relative times of divergence are employed to evaluate the relationship 
between age and diversity among the four major lineages whose distributions broadly overlap on the islands. For 
Hynobius, our analyses are congruent with the previously reported high level of cryptic diversity in morphology and 
allozymes, particularly in species composed of non-sister matrilines. Both species and genetic diversity correlate with 
the relative ages of the lineages. This correlation indicates that the variation in levels of diversity can be explained, to a 
considerable extent, by the hypothesis that older insular lineages have accumulated greater diversity. In addition to the 
Korean Peninsula, H. leechii might have survived in another Pleistocene glacial refugium north of the peninsula and this 
refugium provided a source of colonization after the last glacial maximum. 
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1. Introduction 


An understanding of patterns of biodiversity requires 
the integration of historical and environmental factors 
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(Svenning and Skov, 2005; Hawkins er al., 2006; 
Donoghue, 2008; Jetz and Fine, 2012). Species diversity 
in a region is often considered as the product of net 
diversification rate (speciation rate minus extinction 
rate), time, and dispersal (Ricklefs, 2004; Wiens and 
Donoghue, 2004; Rabosky, 2009). Thus, time is an 
important historical factor when interpreting patterns of 
species richness. The idea that species richness might 
be correlated with how long the constituent lineages 
have been evolving in the area has a long history in the 
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evolutionary time hypothesis (Willis, 1922; Fischer, 
1960; Stebbins, 1974) and the time-for-speciation effect 
(Stephens and Wiens, 2003). These perspectives and their 
variants help explain patterns of species diversity on a 
variety of spatial scales and in various taxonomic groups 
(Borges and Brown, 1999; Hawkins et al., 2005, 2007; 
Pyron and Burbrink, 2009; Roncal et al., 2011), including 
amphibians (Wiens et al., 2006, 2009; Kozak and Wiens, 
2010). 

Lineages of salamanders on Japanese islands exhibit 
quite different levels of species diversity. Five genera 
occur on islands north of a zoogeographical boundary 
known as Watase’s line between the Palearctic and 
Oriental regions (Okada, 1927; Tanaka er al., 1975): the 
hynobiids Salamandrella keyserlingii, Onychodactylus, 
and Hynobius; the salamandrid Cynops pyrrhogaster; 
and the cryptobranchid Andrias japonicus. Among 
them, the distribution of S. keyserlingii is restricted 
to the northeastern margin of the islands, which is 
also the southeastern edge of the current distribution 
of this widespread Eurasian species (Poyarkov and 
Kuzmin, 2008). The other four genera generally 
range widely and overlap extensively on the islands 
(Figure 1). Andrias japonicus is endemic to Japan, 
and it has limited genetic variation throughout its 
range (Matsui er al., 2008). Onychodactylus occurs 
on Japanese islands and the adjacent Asian mainland. 
Of the two species endemic to Japan, O. japonicus 
has a high level of genetic variation and possible 
cryptic taxonomic diversity (Yoshikawa er al., 
2008, 2010a, b, 2012; Poyarkov et al., 2012). Cynops 
pyrrhogaster is endemic to Japan and also has a high 
level of genetic diversity, suggesting the existence of 
cryptic species (Tominaga et al., 2010; Tominaga et al., 
in press). Hynobius occurs on Japanese islands, 
Taiwan Island (China), and mainland Asia. Eighteen 
species of Hynobius are endemic to Japan (Frost, 
2011) and potential cryptic taxonomic diversity 
occurs within several species (e.g., Tominaga et al., 
2005a; Matsui et al., 2006; Nishikawa et al., 2007). The 
occurrence of each of these four genera in Japan likely 
owes to a single colonization event that largely either 
concurred with or after the Miocene separation of the 
islands from mainland Asia (Maruyama et al., 1997; 
Matsui et al., 2008; Yoshikawa et al., 2008; Zhang P. et al., 
2008; Tominaga er al., 2010; Li et al., 2011; Zheng et al., 
2011). 

The various levels of diversity among these lineages 
may be a function of time since colonization. Although 
lineage-ages have been estimated previously, cross- 
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study comparisons are complicated by the employment 
of differing calibration strategies and molecular markers, 
among-lineage rate variation, and innate uncertainty of 
the molecular clock itself (eg. Matsui et al., 2008; Zhang 
P. et al., 2008; Tominaga er al., 2010). The simultaneous 
estimation of ages based on a molecular phylogeny 
containing all these lineages can circumvent this problem 
when using relaxed molecular clock methods. 

Estimations of taxonomic diversity of constituent 
lineages can facilitate evaluations of the relationship 
between age and diversity. Phylogenetic analyses based on 
extensive sampling are available for Japanese lineages of 
Andrias, Onychodactylus, and Cynops (e.g., Matsui et al., 
2008; Yoshikawa et al., 2008; Tominaga et al., in press); 
Hynobius requires evaluation. In Japan, some species 
of Hynobius may contain cryptic taxonomic diversity, 
especially H boulengeri, H. naevius, H. nebulosus, 
and H. yatsui, as evidenced by allozymes, morphology, 
and matrilineal genealogies (e.g., Tominaga et al., 
2005a, b, 2006; Matsui et al., 2006; Nishikawa et al., 
2007). Previous studies involve a relatively small 
taxonomic scale, i.e., a few species. On the islands, great 
morphological similarity among species of Hynobius 
confounds their identification; often locality data are 
necessary (Nishikawa et al., 2007). Moreover, newly 
revealed cryptic species are not always the sister lineages 
of the named species (Donald et al., 2005; Bickford et al., 
2007; Seifert, 2009). A genealogical study involving 
widespread sampling, including not only the Japanese 
species but also species from other areas, may reinforce 
the results of previous studies on Hynobius and provide 
new insights. It will not only facilitate the evaluation of 
diversity of this genus, but also provide an understanding 
of the dispersal history of its members. 

The Korean Peninsula and adjacent areas host six 
species and proposed species of Hynobius (Baek et al., 
2011). Except for H. leechii, all other species occur on the 
southern part of the peninsula only. Within its southern 
range, H. leechii has high levels of mitochondrial DNA 
(mtDNA) diversity (Baek et al., 2011) yet very low 
levels of diversity in allozymes and mtDNA exist in 
the northern range of this species (Zeng and Fu, 2004). 
This “southern richness and northern purity” pattern is 
consistent with the southern Korean Peninsula providing 
a glacial refugium for both animals and plants (Hewitt, 
2000; Kong, 2000; Serizawa et al., 2002; Zhang H. et al., 
2008). Notwithstanding, Zeng and Fu (2004) reported 
a substantial genetic difference between the southern 
(two individuals) and northern forms of H. leechii, thus 
indicating a possible northern refugium outside the 
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Figure 1 Distribution of salamanders on the Japanese islands. Distributions (black) of the species of Hynobius and Cynops pyrrhogaster 
were obtained from IUCN (2011). The ranges (black) of Onychodactylus and Andrias japonicus were mapped following Yoshikawa et al. 


(2008) and Matsui et al. (2008), respectively. 


peninsula (Provan and Bennett, 2008), as proposed for 
rodents (Serizawa et al., 2002; Lee et al., 2008). 

Herein, we first explore the matrilineal history and 
diversity of the genus Hynobius based on comprehensive 
sampling and apply the results to propose taxonomic 
hypotheses. On the basis of this and previous molecular 


analyses of salamanders, we simultaneously estimate 
the relative divergence times to evaluate the relationship 
between lineage-age and current diversity among 
salamanders on the Japanese islands. The species S. 
keyserlingii is not included in this study because its 
restricted insular distribution is considerably isolated from 


No. 4 Yuchi ZHENG et al. 


the ranges of most other lineages. Finally, we examine the 
possibility of a northern glacial refugium for H. leechii 
using increased sample sizes and sampling sites from both 
the northern and southern parts of its distribution. 


2. Material and methods 


2.1 Taxon sampling and DNA sequencing Genealogical 
analyses of Hynobius used a total of 368 individuals from 
30 of the 32 species (Frost, 2011) as the ingroup. Three 
sets of DNA sequence data were available for this genus. 
The first included the entire mitochondrial genome (e.g., 
Zhang et al., 2006; Zheng et al., 2011). Second, data 
from the mitochondrial gene cyt b included 1141 bp, the 
full length (e.g., Matsui et al., 2007; Lai and Lue, 2008; 
Sakamoto er al., 2009; Baek et al., 2011). The third data 
set included two mitochondrial fragments: a /2S—/6S, 
~1075 bp fragment including the intervening tRNA gene 
(GenBank accession No. AY915973-96) and a ~1455 bp 
fragment consisting of the complete ND2 and part of COI 
genes and the tRNA genes between them (AY915925—48) 
(Macey et al., 2005, unpublished data). All these GenBank 
sequences of Hynobius as of February, 2012 were 
included in the analyses. Several individuals with cyt b 
sequence data also had a /2S sequence (Tominaga et al., 
2006) that overlapped substantially with the 7265-165 
fragment; these /2S sequences were included in the 
analysis. We sequenced the mitochondrial genomes of one 
individual each of five species, cyt b from 50 individuals 
representing 23 sampling sites plus three individuals 
without precise locality data, and the /2S—/6S fragment 
of one individual (JQ929919-77). Mitochondrial 
genomes of five species (H. hidamontanus, H. kimurae, 
H. lichenatus, H. nigrescens, and H. tsuensis) were 
sequenced for the first time with the purpose of including 
more distinct full-length sequences in the supermatrix 
approach. The other samples sequenced concentrated 
on mainland species, especially H. leechii. Two species 
of Batrachuperus and Liua were included in the 
outgroup based on the current phylogeny (Peng et al., 
2010; Zheng et al., 2011). Details of the sampling and 
PCR primers were presented in Supplementary Material 
1. No DNA sequences were available from H hirosei and 
H. turkestanicus. 

Sequences of all nine mitochondrial light-strand 
encoded genes were converted into their complementary 
strand. Alignment was conducted with ClustalX v. 1.83 
(Thompson er al., 1997) and checked by eye. Homologies 
of non-coding genes were checked against the secondary 
structures of tRNAs determined using tRNAscan-SE 
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v. 1.21 (Lowe and Eddy, 1997) and the rRNA secondary 
structures of Xenopus laevis (Cannone et al., 2002) and 
the salamander Ambystoma mexicanum (Wuyts et al., 
2004). Amino acid sequences were used to confirm 
homologies for coding regions. Sites of questionable 
homology were excluded from analysis. 


2.2 Molecular phylogenetic analysis of the genus 
Hynobius All 13 protein-coding, two rRNA, and 22 
tRNA mitochondrial genes were analyzed in combination. 
Because no overlaps occurred between sequences of 
some haplotypes, two datasets were analyzed separately, 
each composed of different subsets of taxa based on 
aligned sequences (Supplementary Material 2). In the 
larger dataset, Data-L, all haplotypes had an overlapping 
fragment of cyt b. In the smaller dataset, Data-S, /2S 
overlapped. A six-partition strategy was applied to both 
datasets. A separate partition was defined for each codon 
position from all protein-coding genes, another for each 
rRNA gene, and one partition for the concatenated tRNA 
genes (Mueller et al., 2004; Zhang and Wake, 2009). The 
Bayesian information criterion (BIC) and the corrected 
Akaike information criterion (AICc) implemented in 
jModelTest v. 0.1.1 (Posada, 2008) were used to select 
an evolutionary model that best fit each dataset partition 
(Posada and Buckley, 2004; Tamura et al., 2011). 

Bayesian inference (BI) and maximum likelihood 
(ML) analyses were conducted on both datasets. BI 
was performed using MrBayes v. 3.1.2 (Ronquist 
and Huelsenbeck, 2003) on the CIPRES web portal 
(Miller et al., 2010). Four Monte Carlo Markov chains 
(MCMC) were used to obtain 20 million generations. Six 
independent runs were performed to ensure that analyses 
were not trapped in local optima. Trees were sampled 
every 1000 generations and the last 10 000 sample trees 
were used to construct a majority rule consensus tree and 
the frequency of nodal resolution was termed a Bayesian 
posterior probability (BPP). ML analyses were conducted 
using RAxML v. 7.2.6 (Stamatakis, 2006). The rapid 
hill-climbing algorithm (Stamatakis et al., 2007) was 
used and 200 inferences were executed. Nodal support 
was estimated with nonparametric bootstrap proportions 
(Felsenstein, 1985) involving 1000 replicates. In addition, 
the Shimodaira-Hasegawa (S-H) test (Shimodaira and 
Hasegawa, 1999) was conducted using RAxML to 
evaluate the significance of alternative topologies. 

For H. leechii, an unrooted network approach was 
employed to visualize associations among the haplotypes 
of cyt b. At the 95% confidence level, the statistical 
parsimony analysis was performed using TCS v. 1.21 
(Clement et al., 2000). 
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2.3 Estimating sequence divergence In Hynobius, 
genetic diversity assessments for cyt b and /2S used 
uncorrected nucleotide p-distances calculated in MEGA 
v. 5.05 (Tamura et al., 2011) based on the alignment for 
genealogical analysis. Pairwise distances were calculated 
for all haplotypes and between-group mean distances 
were obtained when necessary. Substantial p-distances 
between H leechii and H. yangi, H boulengeri and H. 
kimurae, and H. amjiensis and H. yiwuensis served as 
conservative interspecific reference points. 


2.4 Molecular clock analysis A total of 38 species 
of salamanders were selected as representatives for 
the clock analysis. Most hypothesized divergence 
events, i.e., nodes in phylogenetic trees, leading to 
Japanese lineages of the genera Andrias, Cynops, 
Hynobius, and Onychodactylus were included; several 
weakly supported events were not included to avoid 
historical uncertainty. In addition to the /2S—/6S and 
cyt b fragments, only the mtDNA fragment ND/—COJ 
(approximately 2040 bp before alignment) was also 
available for all these species as of February, 2012. These 
three fragments were concatenated for analysis. Two 
anuran species, Xenopus tropicalis and Ascaphus truei, 
were selected as outgroup members (Cannatella et al., 
2009; Zhang and Wake, 2009). Details of sampling for the 
molecular clock analysis were provided in Supplementary 
Material 3. We employed the same methods of 
sequence alignment, dataset partitioning, and selection 
of a substitution model as noted above. The sequence 
alignment was presented in Supplementary Material 2. 
We used a reference topology involving 38 sampled 
species of salamanders for the molecular clock analysis. 
This consisted of results that were congruent for 
previous molecular studies and our analyses. Compared 
with our three-fragment dataset, Data-MC, some other 
datasets had more comprehensive taxon sampling for 
particular lineages and/or more loci (1.e., the complete 
mitochondrial genome) (Vieites et al., 2007; Matsui et al., 
2008; Yoshikawa et al., 2008; Zhang P. et al., 2008; Zhang 
and Wake, 2009; Zheng et al., 2011). The likelihood 
ratio test rejected (P < 0.001) the hypothesis that our data 
evolved according to a strict molecular clock. Therefore, 
two widely used, relaxed molecular clock methods 
were adopted: the penalized likelihood approach (PL; 
Sanderson, 2002) and the Bayesian approach developed by 
Thorne, Kishino, and their colleagues (TK; Thorne et al., 
1998; Kishino et al., 2001; Thorne and Kishino, 2002). 
No calibrations were used in the molecular clock 
analysis because we focused on the temporal sequence of 
divergence. Given the poor fossil record for salamanders, 


calibration would have used a large probability 
distribution (Anderson, 2008, 2012; Zhang and Wake, 
2009) and estimated divergence times of this group would 
have varied broadly (e.g., Vieites et al., 2009). Such 
uncertainty would likely obscure the temporal sequence 
of interest. As a result, and also because the uncertainty in 
the root age would be a major source of uncertainty in the 
estimates of node ages (Wiens, 2007), the ingroup root, 
the Cryptobranchoidea—Salamandroidea split, was fixed 
with an arbitrary number. 

For PL, branch lengths in the ML framework 
(Schwartz and Mueller, 2010) were estimated on the 
reference topology using RAxML. The PL analysis was 
performed using r8s v. 1.71 (Sanderson, 2003), in which 
estimates had two decimal places. The ingroup root 
was fixed at the arbitrary value of 500 to facilitate the 
comparison of estimates. The smoothing parameter value 
was set to 16, which was selected from values between 
0.01 and 10 000 by conducting a cross-validation test. To 
assess error levels in estimates, 200 partitioned bootstrap 
replicate datasets were generated (in RAxML) and 
analyzed with the help of Torsten Eriksson’s r8s-bootkit 
available at http://www.bergianska.se/index_forskning ` 
soft.html (Sanderson and Doyle, 2001). In doing this, the 
cross-validation test was performed for each bootstrapped 
dataset. 

The TK analysis followed the manual of Rutschmann 
(2005) and was implemented with PAML v. 4 (Yang, 
2007) and Multidivtime (Thorne and Kishino, 2002). The 
ingroup root was constrained to a value between 4.999 
and 5.001 because calibrations could not be fixed in 
Multidivtime and an ingroup root age prior of between 0.1 
and 10 time units was recommended in the Multidivtime 
readme file. The priors for rate of evolution at the root 
branch and the standard deviation (SD) of the rate were 
set to the same value, which was 1/2 of the mean of ML 
distances between species-pairs descended from the root 
node divided by 5. The priors for the Brownian motion 
constant and the SD of the constant were both set to 1. 
The beta prior on proportional branch depth was set to 
1. The number of samples, cycles between samples, and 
cycles before the first sample were set to values of 10 000, 
100, and 2 000 000, respectively. The analysis was run 
twice and the results were compared to ensure that the 
MCMC reached stationarity. 

As a comparison of proportional time estimates for 
different divergence events, ratios between the estimates 
were calculated. Calculation of the ratio used estimates 
from PL based on Data-MC and mean estimates from the 
TK approach. The 95% confidence interval of the ratio 
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was estimated using either 200 bootstrap replicates (PL) 
or the 10 000 samples of the MCMC (TK). A ratio value 
was first calculated for each bootstrap replicate or sample, 
and then the confidence interval was calculated as either 


bootstrap mean + 1.96 SD (PL; normality not rejected at a 
0.05 level by the Kolmogorov-Smirnov test implemented 
in SPSS v. 12.0) or by sorting the 10 000 resulting values 
and reporting the 250" and 9750" values (TK). 


3. Results 


3.1 Molecular phylogenetic analysis Data-L contained 
263 haplotypes and 15 069 nucleotide sites, of which 
5061 sites were variable, and 3780 were potentially 
parsimony-informative among the ingroup members. 
Among the haplotypes, 23 were based on complete 
or nearly complete mitochondrial genomes and 262 
overlapped cyt b at 637 sites. The exception, Hkim-1, 
contained 370 sites of the overlapping region. Data-S had 
53 haplotypes and 15 069 positions, of which 5081 were 
variable and 3870 were potentially parsimony-informative 
among the ingroup members. A total of 23 haplotypes 
were based on mitochondrial genomes, 46 overlapped by 
2334 positions, and all haplotypes overlapped /2S at 532 
sites. In both datasets the overlapping regions between 
coding genes A7P8 and ATP6 (ten positions), ND4L and 
ND4 (seven positions), and ND5 and ND6 (15 positions) 
were treated as a second codon position. RAxML applied 
one substitution model to all DNA data partitions and the 
GTR+I+G model was used. Models with a proportion 
of invariable sites were mostly selected for individual 
partitions. For BI, a same set of model parameters were 
selected by BIC and AICc for Data-L, but different 
parameters were selected by the two criteria for Data-S. 
Consequently, separate BI analyses of Data-S used 
parameters selected by BIC and AICc and the results were 
compared. Substitution models selected for individual 
dataset partitions were listed in Supplementary Material 4. 
The ML and BI analyses of Data-L produced two very 
similar topologies with most major lineages being well- 
supported. As the only notable difference, H. glacialis 
formed the sister group of H arisanensis + H. sonani 
on the ML tree (Figure 2), while the relationship among 
H. glacialis, H. formosanus, and H arisanensis + H. 
sonani remained unresolved on the BI tree. All the other 
differences involved poorly supported, intraspecific 
relationships, mostly within particular lineages. Both 
approaches resolved two lineages for H. leechii, termed A 
and B. Lineage A was distributed in the southern part of 
the species’ range. In Lineage B, all the three haplotypes 
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from the northern area formed a sublineage (Figure 2, 
Northern Sublineage) that clustered within samples from 
the southern area. 

The ML and BI analyses of Data-S produced nearly 
identical topologies in which most nodes were well- 
supported. Independent BI analyses using model 
parameters selected by BIC and AICc produced identical 
topologies and similar BPPs. Only one difference 
emerged between the ML and BI trees; the ML tree 
(Supplementary Material 5) resolved Hyat-1 as the sister 
group of the lineage containing haplotypes Hhid-G, 
Hbou-100, and Hnae-100. In the BI tree, the relationships 
among Hyat-1, Hnae-1, and the lineage containing 
Hhid-G, Hbou-100, and Hnae-100 were unresolved. 

Weakly-supported alternative topologies were obtained 
for haplotypes Hnae-1, Hneb-G, and Hyat-1 from Data-L 
and Data-S yet these trees were mostly compatible 
(Figure 2; Supplementary Material 5). Thus, some 
samples included only in Data-S were integrated into the 
Data-L tree (Figure 2). Combined, a total of 39 major 
lineages of Hynobius were identified. Some haplotypes of 
H. boulengeri, H. naevius, H. nebulosus, and H. yatsui did 
not cluster with conspecific samples. For H nebulosus, 
the forced unification of all haplotypes from Kyushu 
Island (Hneb-G, Hneb-9, and Hneb-100) in the Data-S 
tree was not rejected by the S-H test (P > 0.05). The same 
result occurred for Hneb-G and Hneb-9 in the tree from 
Data-L (P > 0.05). Our gene trees were similar to those 
of previous studies (Matsui et al., 2007; Lai and Lue, 
2008; Nishikawa et al., 2010; Baek et al., 2011; Li et al., 
2011; Zheng et al., 2011). 

For H. leechii, after excluding sites with missing data, 
the network analysis for cyt b contained 740 sites and 
27 haplotypes. It resolved two groups corresponding to 
Lineage A and Lineage B. In the network corresponding 
to Lineage B (Figure 3), the three northern haplotypes 
clustered together and connected to the interior southern 
haplotypes through 10 unobserved intermediate 
haplotypes. 


3.2 Sequence divergence Pairwise and between-group 
p-distances were estimated for overlapping regions 
of cyt b (637 sites) and /2S (532 sites) (Table 1). The 
interspecific reference distances ranged of 6.31%-12.66% 
(cyt b) and 1.69%—3.01% (12S). Similar levels of 
divergence were observed when Hnae-100, Hyat-1, Hyat- 
2, Hbou-100, Hneb-8, and Group-II (Figure 2) were 
compared with other lineages. These six lineages were 
from the four species that did not cluster into a single 
matriline. 
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Figure 2 The matrilineal genealogy of Hynobius inferred from a maximum likelihood (ML) analysis of the large mitochondrial dataset 
Data-L in which all haplotypes overlapped the fragment for cyt b. Nodes with ML bootstrap proportions (BP) > 90 and Bayesian posterior 
probabilities (BPP) > 95 are indicated as asterisks. Nodes with 90 > BP > 70 and BPP = 95 and nodes with 90 > BP > 70 and 95 > BPP > 80 
are marked by closed and open circles, respectively. Vertical bars indicate species designation or lineage assignment. Dashed branches 
indicate mapped samples from small dataset Data-S where the fragments of 72S overlapped. Arrow indicates a node not recovered in the 


Bayesian inference analysis. Haplotype names correspond with those in Supplementary Material 1. 
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3.3 Proportional time estimates Proportional time 
estimates for four splitting events were compared, 
including the basal split in Hynobius, the split between 
the Japanese lineage of Onychodactylus and O. zhaoermii 
from northeastern China, between C. pyrrhogaster 
and C. ensicauda, and between A. japonicus and A. 
davidianus. For PL, branch lengths were estimated with 
the GTR+I+G model (Supplementary Material 4). PL and 
TK approaches yielded similar estimates for these events 
and the same three-level temporal sequence of the events 
(Figure 4). The basal split in Hynobius predated all others. 
Splits within Onychodactylus and Cynops occurred near 
simultaneously and were relatively recent events, ranging 
from 0.54 to 0.67 times that of Hynobius. The split for 
Andrias ranged from 0.22 to 0.24 of that for Hynobius, 
and thus it was the youngest. Differences between these 
divergences time were statistically significant; the 95% 
confidence interval of a ratio between two levels of time 
was less than one. Several divergences between Japanese 
species of Hynobius were temporally close to the basal 
split of this genus in the PL and TK analyses (Figure 4). 


4. Discussion 


4.1 Age-diversity relationship in salamanders of the 
Japanese islands Time is critical for understanding the 
drivers of diversity (Wiens er al., 2009; Rabosky, 2012; 
Kozak and Wiens, 2012). Among the four lineages of 


Hlee-18 


Lineage B 


Table 1 Estimated nucleotide distances among salamanders of the 
genus Hynobius. Haplotype and lineage names correspond with 
names in Figure 2 and Supplementary Material 1. 


Nucleotide p-distance (%) 


Comparison 
cyt b (637 sites) 12S (532 sites) 

H. leechii — H. yangi 6.31° 1.69° 
H. boulengeri — H. kimurae 8.95° 3.01° 
H. amjiensis — H. yiwuensis 12.66° 1.72° 
Hnae-1 — other haplotypes > 10.20 = 1.88 
Hnae-100 — other haplotypes — > 2.44 
Hyat-1 — other haplotypes > 9.89 > 1.88 
Hyat-2 — other haplotypes > 8.32 = 2.63 
Hbou-100 — other haplotypes = > 2.26 
Hneb-8 — H. tokyoensis 13.44 — 
Hneb-8 — Group-II ODER — 
Group-l — H. tokyoensis 11.19" — 
Hneb-G — Hneb9 7.85 — 
Hneb-G — (Hneb9 + Hneb100) — 1.22° 


* Between-group mean distances. 


Japanese salamanders studied, lineage-age estimates are 
positively related to species richness and phylogenetic 
diversity. 

In East Asia, the ancestral distribution of extant 
Hynobius is estimated to occur on land masses of the 
current Japanese islands (Li et al., 2011). Thus, it is likely 


Figure 3 Network for cyt b haplotypes from Lineage D of Hynobius leechii and sampling sites of this species. Lines represent single 
mutational changes. Ovals are sampled haplotypes. Closed circles represent haplotypes that are necessary intermediates but not encountered. 
Shaded ovals are haplotypes found in the northern sampling sites indicated by triangles. Occurrences of haplotype Hlee-92 are indicated by 
closed triangles. Haplotype names correspond with those in Figure 2 and Supplementary Material 1. Two haplotypes used in the genealogical 
analysis are identical after excluding positions with missing data. Locality data for three specimens are not precise enough to be included in 


the map (Supplementary Material 1). 
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that this area was colonized by Hynobius no later than the 
basal split (node H, Figure 4) among East Asian members 
of the genus, which is dated at about 11—23 million years 
ago (Zheng et al., 2011). As the oldest lineage, this genus 
gradually diverged long ago to form 18 species plus some 
possible cryptic species in Japan (Figures 2, 4). In contrast 
to Hynobius, the genera Andrias, Onychodactylus, and 
Cynops each have one endemic lineage (one or two species) 
on the islands. All of them nest within relatives occurring 
outside the islands (Matsui et al., 2008; Yoshikawa et al., 
2008; Zhang P. et al., 2008; Tominaga et al., 2010). 
Consequently, the divergence times between these 
endemic lineages and their congeneric sister taxa (C and 
A, Figure 4) or close relatives (O, Figure 4; see below) 


TK 


AIH = 0.24 (0.16-0.35) 
O/H = 0.58 (0.43-0.79) 
C/H = 0.67 (0.48-0.91) 
AIO = 0.41 (0.27-0.62) 
AIC = 0.36 (0.24-0.55) 


PL 


AIH = 0.22 (0.16-0.29) 
O/H = 0.58 (0.42-0.74) 
C/H = 0.54 (0.41-0.73) 
AIO = 0.38 (0.27-0.52) 
AIC = 0.40 (0.27-0.53) 
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are used here as indicators of how long each lineage has 
been evolving on land masses of the islands. 

The lineages of Onychodactylus and C. pyrrhogaster 
have moderate periods of evolutionary time and 
corresponding levels of diversity. Date estimates for 
splitting events O and C are significantly younger in being 
0.54 to 0.67 times less than that for H (Figure 4). In O. 
japonicus sensu lato, using complete sequences of cyt b, 
Yoshikawa et al. (2008) recognized four considerably 
differentiated lineages (p-distances 5.5%—9.6%). They 
suggested the presence of cryptic species. In support 
of this, extensive genetic and morphological variation 
occurs in this species and several candidates of cryptic 
species and one new species, O. nipponoborealis, were 


Cryptobranchus alleganiensis 
A Andrias davidianus 
Andrias japonicus 
O Onychodactylus zhaoermii 
Onychodactylus japonicus sensu lato 
Paradactylodon mustersi 
Pachyhynobius shangchengensis 
Liua tsinpaensis 
Batrachuperus londongensis 
Hynobius kimurae 
Hynobius arisanensis 
Hynobius nigrescens 
Hynobius lichenatus 
Hynobius tokyoensis 
Hynobius quelpaertensis 
Hynobius leechii 
Hynobius yangi 
Hynobius nebulosus 
H Hynobius tsuensis 
Hynobius hidamontanus 
Hynobius yiwuensis 
Hynobius amjiensis 
Hynobius guabangshanensis 
Hynobius chinensis 
Hynobius maoershanensis 
Hynobius retardatus 
Eurycea bislineata 
Ambystoma mexicanum 
Salamandrina terdigitata 
Lyciasalamandra flavimembris 
Echinotriton andersoni 
Taricha rivularis 
Euproctus platycephalus 
Triturus cristatus 
Cynops orientalis 
Cynops orphicus 
Cynops pyrrhogaster 
C Cynops ensicauda 


Figure 4 The time-calibrated tree of salamanders based on the mitochondrial dataset Data-MC containing 3693 nucleotide positions. 
Proportional times for nodes were estimated using the Bayesian approach developed by Thorne, Kishino, and their colleagues (TK). Gray 
bars through the nodes indicate 95% highest posterior densities for estimates. Numbers in parentheses are 95% confidence intervals. PL = the 
penalized likelihood approach. Bold species names indicate species that occur on the Japanese islands north of Watase’s line. The outgroup is 


not shown. 


298 Asian Herpetological Research 


Vol. 3 


reported in follow-up studies (Yoshikawa et al., 2010a, b, 
2012; Poyarkov et al., 2012). The same magnitude of 
difference for the complete cyt b sequences also occurs 
in C. pyrrhogaster, based on limited samples of this 
wide-ranging species (Tominaga et al., 2010). Based on 
a comprehensive sampling, Tominaga et al. (in press) 
recently suggested that they were actually four species. 
Andrias japonicus has the shortest period of evolutionary 
time and lowest level of genetic diversity. Date estimates 
for the split of Andrias (A) range from 0.22 to 0.24 times 
of the estimate for Hynobius, and this is significantly later 
than the splits of Onychodactylus and Cynops. Compared 
with its sister species in mainland Asia, A. japonicus 
exhibits less genetic divergence (Murphy et al., 2000; 
Matsui et al., 2008). These correlations indicate that the 
varying levels of diversity among lineages of salamanders 
on the Japanese islands likely owe, to a considerable 
extent, to the amount of time that each lineage has been 
evolving in the region. 

Our molecular clock analyses cover a time-span of 
about 200-300 million years (Roelants er al., 2007; 
Vieites et al., 2007, 2009). Substitutions between some 
highly divergent sequences are likely severely saturated 
(ep, Zheng er al., 2011). However, this does not imply 
that the evolutionary time-order revealed from the 
analysis is unreliable. In the molecular clock approach, 
divergence times between organisms relate to numbers 
of substitutions accumulated in sequences. Saturation 
due to multiple hits usually results in greater percentage 
of underestimated substitutions between more divergent 
sequences (Nei and Kumar, 2000; Arbogast et al., 2002). 
As proportional time estimates were compared as ratios in 
this study, the effect of saturation would likely compress 
the difference between divergence dates, rather than result 
in a Statistically significant order of time. 

Poyarkov et al. (2012) recently described a new 
species, Onychodactylus zhangyapingi, which may 
complicate our analysis. Onychodactylus zhangyapingi 
occurs in northeastern China and is the sister group 
of Japanese congeners (Poyarkov et al., 2012). This 
relationship suggests that we likely overestimated the 
relative age of Japanese Onychodactylus (O). The reason 
of the exclusion of this species in our final analysis is 
that the available data for this new species are far less 
than those in dataset Data-MC. Nevertheless, the splitting 
event O between O. zhaoermii and (O. zhangyapingi 
+ Japanese Onychodactylus) and the split between O. 
zhangyapingi and Japanese Onychodactylus are close to 
each other and both are relatively deep (Poyarkov et al., 
2012). Therefore, our overestimate is likely insignificant 


to the age-diversity pattern (Figure 4). 


4.2 Cryptic diversity in Japanese Hynobius In 
addition to lineages corresponding to the 17 (of 18) 
species sampled, six other major lineages were identified 
including Hnae-100, Hyat-1, Hyat-2, Hbou-100, Hneb-8, 
and Group-II (Figure 2; Supplementary Material 5). The 
genetic distances between them and other lineages were at 
least similar to the reference interspecific distances of the 
genus (Table 1). The two major lineages of H naevius, 
Hnae-1 and Hnae-100, were not sister lineages and they 
corresponded to lineages identified by Tominaga et al. 
(2006). Lineage 1 was from northwestern Kyushu Island 
(Hnae-1) and Lineage 2 from the other populations 
(Hnae-100). Because the type locality of H. naevius 
was estimated to be in the northwestern Kyushu Island 
(Tominaga and Matsui, 2007), the other populations 
may be a cryptic species. Consistent with this finding, 
these two groups differ by two of 20 allozyme loci (fixed 
alleles at s/DH-A and nearly fixed at m4COH-A) and by 
two-dimensional scaling of allozyme data based on 19 
(of 20) polymorphic loci (Tominaga et al., 2005a). Two- 
dimensional scaling analyses of morphology also separate 
the lineages (Tominaga et al., 2005b; Tominaga and 
Matsui, 2008). 

Three major lineages occur among samples of H. 
yatsui: the lineage from Kyushu Island, Hyat-2, and 
Hyat-1, of which the former two are sister-lineages 
(Figure 2; Supplementary Material 5). The former lineage 
contains all haplotypes on the island and this includes 
the type locality of H. yatsui (Tominaga and Matsui, 
2008). This lineage, named KYUSHU or Kyushu-B, is 
clearly distinguishable from non-Kyushu populations by 
allozymes (e.g., loci PGDH-A and PGM-C), morphology, 
and matrilineal studies (Tominaga et al., 2005a, b, 2006; 
Tominaga and Matsui, 2008). The non-Kyushu Island 
morphological groups and genetic lineages, including 
Hyat-1 and Hyat-2, suggest cryptic diversity and this 
requires further study. 

Samples of H. boulengeri comprise two distantly 
related lineages. The lineage represented by haplotypes 
Hbou-! and Hbou-101 corresponds to the true H. 
boulengeri, as restricted to Honshu Island by Nishikawa 
et al. (2007) based on morphological and allozymic 
variation. The other lineage, Hbou-100, may represent 
either an unnamed species (Nishikawa et al. 2007) or 
H. hirosei, which was removed from the synonymy of 
H. boulengeri by Nishikawa et al. (2007). The precise 
collecting locality of Hbou-100 remains unknown and 
this precludes a determination. 

Hynobius nebulosus has three haplotypes from 
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near the type locality of Nagasaki: Hneb-G, Hneb- 
9, and Hneb-100. Other, distantly related samples 
form two major lineages that are not sister groups: 
Group-II and Hneb-8 (Figure 2). Lineages represented 
by the three haplotypes may be sister-groups (S-H 
test, P > 0.05) and the genetic distances between 
them are less than most values for the reference 
species (Table 1). Therefore, we follow Matsui er al. 
(2006) in assigning all three haplotypes to H. nebulosus. 
Allozyme analyses for H. nebulosus sensu lato identify 
three candidate cryptic species: eastern, montane, and 
Chugoku groups (Matsui er al., 2006). Our Group-II 
corresponds to the eastern group. Because the sampling 
site (Gobo-shi) of lineage Hneb-8 was not included in the 
study of Matsui et al. (2006), we cannot determine if this 
is another of their potential cryptic species. 

Underestimations of variation may obscure our 
understanding of evolutionary processes (Olsson et al., 
2005; Kaliontzopoulou et al., 2011). Our analyses 
reinforce the high level of cryptic diversity within 
Japanese Hynobius as revealed in previous morphological 
and allozymic studies, particularly by the recognition 
of non-monogenealogical (sensu Murphy and Méndez 
de la Cruz, 2010; Gao et al., 2012) species. This finding 
supports the importance of having dense taxon sampling 
in genealogical and taxonomic studies (eg. Olsson et al., 
2005). A cryptic species is not necessarily the sister group 
of the species from which it differs (Donald et al., 2005; 
Bickford et al., 2007; Seifert, 2009), especially for groups 
such as amphibians that exhibit conserved morphological 
evolution and extensive homoplasy (Parra-Olea and 
Wake, 2001; Fouquet er al., 2007). 


4.3 Northern glacial refugium of Hynobius leechii 
Multiple refugia might have served H /eechii during 
Pleistocene glacial cycling. This species currently occurs 
on the Korean Peninsula and adjacent northeastern China 
with a general north-south distribution. Substantial genetic 
diversity occurs on the southern peninsula yet with a 
similar sampling area to the southern one, extremely 
low diversity occurs northwards. The genealogy reveals 
a southern origin of the northern form (Figure 2). This 
pattern indicates that the southern peninsula provided a 
glacial refugium for H. /eechii (Hewitt, 2000; Provan and 
Bennett, 2008), and this has been indicated for a variety 
of organisms (e.g., Kong, 2000; Serizawa er al., 2002; 
Zhang H. et al., 2008). 

For H. leechii, another refugium probably existed north 
of the peninsula. Several observations are consistent with 
this notion. First, among the three northern haplotypes, 
Hlee-92, the only haplotype found in most sampling sites 
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(Supplementary Material 1), is widespread throughout 
the northern sampling area (Figure 3). The two other 
haplotypes, Hlee-G1 and Hlee-91, differ from the former 
by one or two nucleotide substitutions only (Figure 3). 
Their distributions are far more restricted geographically 
and this pattern infers a local origin. These distributions 
suggest a postglacial range expansion (Provan and 
Bennett, 2008). Second, on the network and genealogy, 
these haplotypes substantially differ from the most closely 
related interior southern haplotypes, in the case of the 
network by at least ten unobserved haplotypes involving 
11 substitutions. A p-distance of 1.49% based on the 740 
sites of cyt b separates them. It is unlikely that this level 
of divergence accumulated during a northwards post- 
last glacial maximum (LGM) dispersal of this species 
from the Korean Peninsula refugium (Canestrelli et al., 
2006; Provan and Bennett, 2008). Given the broadly 
used evolutionary rate of 0.64%-—1.00% divergence per 
million years per lineage for vertebrate cyt b gene (e.g., 
Tominaga et al., 2010), these haplotypes should have 
diverged long before the LGM, which occurred about 
0.02 Ma. The northern form probably persisted in situ 
throughout the LGM. These results agree with the notion 
that postglacial colonization of high latitude regions can 
be from local sources within or close to the regions, rather 
than from the more distant, southern refugium (Stewart 
and Lister, 2001; Pearson, 2006). 
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